/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/const.hxx"
#include "uhm/util.hxx"

#include "uhm/object.hxx"

#include "uhm/operation/scheduler.hxx"
#include "uhm/operation/mesh.hxx"
#include "uhm/operation/element.hxx"

#include "uhm/mesh/node.hxx"
#include "uhm/mesh/element.hxx"

#include "uhm/matrix/uhm/matrix.hxx"

#include "uhm/mesh/mesh.hxx"

extern "C" {
#include "defs.h"
#include "struct.h"
#include "macros.h"
#include "rename.h"

  extern void SetUpGraph(GraphType *, int, int, int, 
			 idxtype *, idxtype *, idxtype *, idxtype *, int);
  extern int idxsum(int, idxtype *);
  extern void InitRandom(int);
  extern void AllocateWorkSpace(CtrlType *, GraphType *, int);
  extern void FreeWorkSpace(CtrlType *, GraphType *);
  extern void MlevelEdgeBisection(CtrlType *, GraphType *, int *, float);
  extern void SplitGraphPart(CtrlType *, 
			     GraphType *, GraphType *, GraphType *);
  extern void GKfree(void **,...); 
}

//#define UHM_BUILD_TREE_MEASURE_TIME 1


namespace uhm {
  static bool build_tree_var_4_internal(std::vector< Element > *orphan,
					Mesh m, Element parent, int is_parallel,
					CtrlType *ctrl, GraphType *graph);
  
  // --------------------------------------------------------------
  // ** coarsening 
  bool build_tree_var_4(Mesh m) {

#ifdef UHM_BUILD_TREE_MEASURE_TIME 
    printf(":: begin :: build_tree_var_4\n");
    double t, time;
#endif

    assert(m && mesh_valid(m));
    { // ** multi level coarsening mesh
      // Start :: collect orphans
      std::vector< Element > *orphan = new std::vector< Element >;
      { // ** preprocess                                              
	// After initial processing of nodes, leaves have separated nodes
#ifdef UHM_BUILD_TREE_MEASURE_TIME 
	t = timer();
#endif

	Scheduler_ s;
	s.load(m);

#pragma omp parallel 
        {
#pragma omp single nowait
          s.execute_elements_par(&(op_update_connectivity), true);
        }

#ifdef UHM_BUILD_TREE_MEASURE_TIME 
	time = timer() - t;
	printf(":: preprocess on mesh %lf\n", time);
#endif
      }

      {
#ifdef UHM_BUILD_TREE_MEASURE_TIME 
	t = timer();
#endif

	std::map< int, Element_ >::iterator mit;
	for (mit=m->elements.begin();mit!=m->elements.end();mit++) 
	  if (mit->second.is_orphan()) 
	    orphan->push_back(&mit->second);
        

#ifdef UHM_BUILD_TREE_MEASURE_TIME 
	time = timer() - t;
	printf(":: adding orphan in mesh %lf\n", time);
#endif
      }
      // ---------------------------------------------------------------
      // *** bisection of orphan
      if (orphan->size() > 1) {
	Element parent = m->add_element();
	int n = orphan->size(), wgtflag = 1;

	std::vector< int > xadj, adjncy, adjwgt;
#ifdef UHM_BUILD_TREE_MEASURE_TIME
        t = timer();
#endif

        build_graph(orphan, parent,
                    &xadj, &adjncy, &adjwgt);
        
#ifdef UHM_BUILD_TREE_MEASURE_TIME
        time = timer() - t;
        printf(":: building graph of the mesh %lf\n", time);
#endif
        
	int tvwgt, tpwgts2[2], *label, *where;
	float ubfactor;
	CtrlType ctrl;
	GraphType graph, lgraph, rgraph;
	
	// create graph object to interface to METIS
#ifdef UHM_BUILD_TREE_MEASURE_TIME
        t = timer();
#endif
        
        SetUpGraph(&graph, OP_OEMETIS, n, 1, &xadj[0], &adjncy[0],
                   NULL, &adjwgt[0], wgtflag);
        
#ifdef UHM_BUILD_TREE_MEASURE_TIME
        time = timer() - t;
        printf(":: create metis graph %lf\n", time);
#endif

	// use default control parameter                
	ctrl.CType   = PMETIS_CTYPE;
	ctrl.IType   = PMETIS_ITYPE;
	ctrl.RType   = PMETIS_RTYPE;
	ctrl.dbglvl  = PMETIS_DBGLVL;
	
	ctrl.optype    = OP_PMETIS;
	ctrl.CoarsenTo = 20;
	ctrl.maxvwgt   = 1.5*(idxsum(n, graph.vwgt)/ctrl.CoarsenTo);
 
	InitRandom(-1);

#ifdef UHM_BUILD_TREE_MEASURE_TIME
        t = timer();
#endif

#pragma omp parallel
	{
#pragma omp single nowait
	  {
	    AllocateWorkSpace(&ctrl, &graph, 2);
	    build_tree_var_4_internal(orphan, m, parent, true, 
				      &ctrl, &graph);
	    FreeWorkSpace(&ctrl, &graph);	
	  }
	}

#ifdef UHM_BUILD_TREE_MEASURE_TIME
        time = timer() - t;
        printf(":: recursive partitioning %lf\n", time);
#endif
        
      }
      delete orphan;
    }

    { // ** update generation
#ifdef UHM_BUILD_TREE_MEASURE_TIME
      t = timer();
#endif
      Scheduler_ s;
      s.load(m);
      //      s.execute_leaves_seq(&(op_restore_connectivity));
      //      s.execute_elements_par(&(op_update_connectivity), true);
      s.execute_tree(&op_update_generation, true);

#ifdef UHM_BUILD_TREE_MEASURE_TIME
      time = timer() - t;
      printf(":: post process of mesh %lf\n", time);
#endif

    }
#ifdef UHM_BUILD_TREE_MEASURE_TIME
    printf(":: end :: build_tree_var_4\n");
#endif

    return true;
  }
  
  bool build_tree_var_4_internal(std::vector< Element > * orphan, 
				 Mesh m, Element parent, int is_parallel,
				 CtrlType *ctrl, GraphType *graph) {

    int n=graph->nvtxs, tvwgt, tpwgts2[2];
    CtrlType  ctrl_elt[2];
    GraphType graph_elt[2];
    
    tvwgt = idxsum(n, graph->vwgt);
    tpwgts2[0] = tvwgt/2;
    tpwgts2[1] = tvwgt-tpwgts2[0];

    MlevelEdgeBisection(ctrl, graph, tpwgts2, 1.0);
    SplitGraphPart(ctrl, graph, &graph_elt[0], &graph_elt[1]);

    GKfree((void**)(&graph->gdata), 
	   (void**)(&graph->rdata), 
	   (void**)(&graph->label), LTERM);

    for (int i=0;i<2;++i) {

      if (graph_elt[i].nvtxs > 3) {
	Element elt = m->add_element();
	
	elt->set_parent( parent );
	parent->add_child( elt );

	ctrl_elt[i] = *ctrl;

	if (graph_elt[i].nvtxs > UHM_DISSECTION_TASK_SIZE && 
	    is_parallel) {
	  
#pragma omp task firstprivate(i, elt, ctrl_elt, graph_elt)
	  {
	    AllocateWorkSpace(&ctrl_elt[i], &graph_elt[i], 2);
	    build_tree_var_4_internal(orphan, m, elt, true,
				      &ctrl_elt[i], &graph_elt[i]);
	    FreeWorkSpace(&ctrl_elt[i], &graph_elt[i]);
	  }
	} else {
	  build_tree_var_4_internal(orphan, m, elt, false,
				    ctrl, &graph_elt[i]);
	}

      } else {
	if ( graph_elt[i].nvtxs ) {
	  Element elt;
	  
	  switch ( graph_elt[i].nvtxs ) {
	  case 1: 
	    elt = parent;
	    break;
	  default:
	    elt = m->add_element();
	    elt->set_parent( parent );
	    parent->add_child( elt );
	    break;
	  }
	  
	  for (int j=0;j<graph_elt[i].nvtxs;++j) {
	    orphan->at( graph_elt[i].label[j] )->set_parent( elt );
	    elt->add_child( orphan->at( graph_elt[i].label[j] ) );
	  }


	  GKfree((void**)(&graph_elt[i].gdata), 
		 (void**)(&graph_elt[i].rdata), 
		 (void**)(&graph_elt[i].label), LTERM);
	}
      }
    }
    /*
    if (is_parallel) {

#pragma omp taskwait

    }
    */
    return true;
  }
}
